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Current-voltage characteristics of diluted Josephson-junction arrays: scaling behavior 

at current and percolation threshold. 
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Dynamical simulations and scaling arguments are used to study the current-voltage (IV) charac- 
teristics of a two-dimensional model of resistively shunted Josephson-junction arrays in presence of 
percolative disorder, at zero external field. Two different limits of the Josephson-coupling concen- 
tration p are considered, where p c is the percolation threshold. For p > p c and zero temperature, 
O ' the IV curves show power-law behavior above a disorder dependent critical current. The power- 

law behavior and critical exponents are consistent with a simple scaling analysis. At p c and finite 
f"S !' temperature T, the results show the scaling behavior of a T = superconducting transition. The re- 

sistance is linear but vanishes for decreasing T with an apparent exponential behavior. Crossover to 
non-linearity appears at currents proportional to T 1+l>T , with a thermal-correlation length exponent 
1 vt consistent with the corresponding value for the diluted XY model at p c . 
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I. INTRODUCTION 



There has recently been an interest, both experimentali-i and theoreticali~0, in the resistive behavior of diluted 
Josephson-junction arrays .-(JJA). These systems provide a useful model for several universal transport properties of 
£SJ ■ granular high-T c materialstHI and can also be physically realized as superconducting arraysEm or wire networks^ with 
t-H I accurately controlled parameters using microfabrication techniques. Although most investigations have considered 
the combined effects of disorder and magnetic fields which leads to the vortex glass stateo there are also interesting 
o> ■ questions even in the absence of an external field as a result of disorder. The effects of percolative disorder on the 
resistive transition andjionlinear current-voltage (IV) relation, at finite temperatures, have been studied in artificial 
two-dimensional arraysQ and also in the context of the high-T c oxidealj. For an ordered two-dimensional JJA, which 
is isomorphic to the XY model, it is well-known that the resistive transition is in the Kosterlitz-Thouless (KT) 
P3 ■ universality class where the low-temperature phase has a nonlinear IV relation, V cx I a ^ T \ due to current induced 
vortex-pair unbinding. The temperature dependent exponent a(T) decreases with increasing temperature and reaches 
1 1 a(T c ) — 3 at the transitiontS. In presence of percolative disorder but well above the percolation threshold, this 
behavior still persists with a broaden resistive transition and disorder renormalized critical temperature T c (p) and 
exponent a(T,p), where p is the fraction of superconducting grainsa. Right at the percolation threshold p c , however, 
as for the XY modelEZI, the transition temperature is expected to vanish but a divergent correlation length £ ~ y-^T 
and thermal critical exponent vt can still be defined for T > 0. A natural question arises as to what should be 
the behavior of the IV curves for T > in this limit since the correlation length is finite. In fact, one expects the 
ion length as T — > to have important consequences at finite T for the nonlinear behavior of 
The linear resistance is expected to be finite for any T > but crossover to non-linearity can 
appear for increasing currents due to the finite correlation length with a universal behavior determined by the T = 
transition. It is important to understand this behavior in detail because, in two dimensions and in presence of an 
external magnetic field, an additional T — (vortex-glass like) transition is expectedu near p c which should also have 
a similar effect but it could be difficult to separate these two effects experimentally. Also, in principle, the universal 
behavior at p c may be useful to identify the proximity to the percolation threshold in systems where p c is uncertain 
as in granular materials modeled as diluted JJA. 

There are also other interesting effects of disorder on the current- voltage relation which appear in a different limit 
and are not well understood. Above the percolation threshold p c and zero temperature, the IV characteristics shows a 
power-law behavior, V ~ (J — Jc) a , above a critical current density J c . For example, an exponent a ~ 3 was obtained 
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from numerical simulations by Xia and LeathtJ in d = 2. Clearly, this behavior at T — and finite current density J c 
is unrelated to the condition a(T c ) = 3 at the KT transition mentioned above and requires different considerations. 
In a linearized model of the critical current as a phase-coherence breakdown phenomenonQ, it has been argued that 
the critical current of a diluted array vanishes in the thermodynamic limit. In this case, the apparent finite J c and 
the exponent a, presumably universal, are in fact an artifact of finite system sizes. However, it is possible that 
nonlinearities in the array problem could lead to a finite value even in the large system limitLru. In the latter scenario 
J c is finite but is expected to vanish at the percolation threshold as J c ~ (p—p c ) v with an exponent v — Vp(d— 1), where 
v p is the percolation exponent and d is the dimensionality of the system. Unlike the behavior at finite temperatures, 
little is known however about the precise value of a and the universality class of the current-induced transition leading 
to the power-law behavior for p > p c , under the assumption that there is in fact a well defined J c . Nevertheless, it is 
important to understand these features in detail in order to be able to separate these effects, which arise only from 
disorder, from others induced by a combination of disorder and strong magnetic fields as in the case of the vortex 
glass which also leads to values of a of the same orderEj at T — 0. Indeed, in experimental systems, these effects can 
appear simultaneously and it is possible that some of the power-law behavior observed in the IV characteristics of 
granular high-T c materials that tend to be regarded as a manifestation of vortex glass behavior is actually a result of 
intrinsic geometrical disorder and would persist even in the absence of the external field. 

In this work, we present dynamical simulations of the IV characteristics of diluted two-dimensional JJA, at zero 
external field, in the two different limits discussed above which allow for a scaling analysis in terms of a single length 
scale, (i) At T — and p > p c , we study the power-law behavior in the IV curves above J c as resulting from a current 
driven dynamical transition. We discuss scaling arguments for this behavior which leads to a = vj(z + 1), where vj 
is a correlation length exponent, z a dynamical exponent and compare with numerical results, (ii) At the percolation 
threshold p = p c and finite T we study the scaling of the IV curves as resulting from a T = superconducting 
transition. The results are consistent with a linear resistance at finite temperatures and nonlinear behavior appearing 
at current densities proportional to T l+UT , where iAp is the thermal-correlation length exponent for the diluted XY 
model at p c . 



II. MODEL AND SIMULATION 



We consider a model of resistively shunted JJA consisting of coupled superconducting islands located at the nodes 
of a square network with Josephson and normal currents flowing between them. The nodes. BxeJocated at r — mx + ny 
with unit lattice constant. The current 7 M (r) flowing between r and r + fi, is modeled aa3 - E£lEJ 

7 M (r)=7? jS in(A^(r,i)) + + ^(r,i). (1) 

Here A^0(r, t) is a discrete gradient of the superconducting phases £?(r, t), 7° is the critical current of the junctions, 
and i? rM is a shunt resistance between the islands. The white noise random variable ^(r, t) represents thermal 
Johnson fluctuations in the current with covariance, 

2kT 

(77 M (r, t) v (r', t')) = —5 r , r ,6^,6(t - t'). (2) 

We assume that disorder affects only the Josephson coupling. In the absence of the coupling, 7° = 0, and so only 
normal current flow between the neighboring grains. Disorder effects in R T ^ aie-,assumed to be less important. This 
model has been applied to-the study of transport properties of hiahpTc oxideaia and, in its site dilution version, to 
composite superconductors]. For proximity coupled artificial arraysoB, this model is also a reasonable approximation 
since, as pointed out in Ref. ^, the underlying normal-conducting layer over which the superconducting grains are 
deposited provides a roughly constant R^. Dilution of the junctions is introduced by taking critical currents 7° = 
with probability p or 7^ = with probability 1 — p and constant shunt resistance, R^ — Rq. After combining Eqs. (^) 
with current conservation at each node, A M • 7 A1 (r) = I ext (r), the dynamical equation of motion for 8(r,t) becomes, 

0(r, *) = -£ G(r, r') (r') - A, ■ [7 r % sin(A^(r', t)) + ^(r', t)) } , (3) 

r' 

with G(r, r') the d = 2 lattice Green function. Dimensionless quantities are used with time in units of tj = h/2eR I , 
currents in units of 7o, voltages in units of i?o7o and temperatures in units of hIo/2ekB- We choose periodic boundary 
conditions along the x-direction and open boundary conditions along the y-direction. The array has 7x7 bonds, 
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corresponding to L x (L + 1) nodes. The total current / is injected uniformly along the y-direction with i ext (m, n) = 
J(8n,i ~ fin,L+i) where J is the current density J = I/L . Eqs--{g) are solved with the second order Runge-Kutta- 
Helfand-Greenside algorithm for stochastic differential equationsBa with a time step of At = O.lr,/. Temporal averages 
are taken over a time of 2000tj after a transient time of 500tj. The matrix multiplication by G(r,r') is performed 
by means of a fast Fourier transform and cyclic reduction algorithm as used in Ref. The voltage drop along 

the y-direction is given by 

1 L 

V=-^2(e(m,n = L + l)-6(m,n = l)) (4) 

m— 1 

in dimensionless units, where (. . .} is a time average, and the average electric field is given by E = V/L . 

III. RESULTS AND DISCUSSION 

Figure 1 shows some of the IV characteristics obtained numerically at T = 0. Our averaging time is a factor of two 
larger than the one used in Ref. || to study a site dilution version of the same model in d = 3. As a test of the numerical 
method we included in Figure 1 the calculation for p = 1 where the array behaves as a single junction with a critical 
current density J c — 1. Above the percolation threshold p c — 1/2, there is an infinite cluster of superconducting 
junctions through the system and an apparent finite critical current density J c below which the voltage is very small. 
For p < pc, only isolated finite clusters occurs and the resistance is Ohmic for small currents. 
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FIG. 1. Current-voltage characteristics for a L — 64 array, at zero temperature, as a function of the Josephson-coupling 
concentration p 

Figure 2a shows the behavior of the critical current density J c and linear resistance Rl = limj^o E/J as a function 
of p at zero temperature, for the largest system size used in the simulations, L = 64. J c decreases with p and 
presumably vanishes at the percolation threshold p c while Rl is nonzero only for p < p c and also appears to jwanish 
at p c . Unfortunately, closer to p c our data is not accurate enough to test the expected power law behavioioEl for 
the critical current J c {p) — c {p — p c ) v ', with v = v p (d — 1) = 4/3, and similar behavior for the linear resistance 
Rl ~ (jpc — p) s , with s ~ 1.3. 

The exponent a obtained from a power-law fit to the current-voltage characteristics just above J c , E ~ (J — Jc)^\ 
as a function of p is indicated in Fig. 2b. At p = 1 the system behaves as a single junction and a — 1/2 exactlyEil 
which agrees with the numerical simulation. In presence of disorder for decreasing p, this exponent jumps to a roughly 
constant value, a ~ 2.5(2), for the largest system size L = 64. For the smaller system size its p-dependence is more 
significant but we consider this as an effect of small system sizes where the true asymptotic limit has not been reached. 
The value of a for L = 32 is in fair agreement with other simulationsa at a fixed value of p = 0.9 for comparable 
system sizes even though the model used differs from ours in that dilution affects both the Josephson coupling and 
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the associated shunt resistance simultaneously. We believe that the difference between the models should not affect 
the behavior of a far above p c where there is an infinite cluster of superconducting junctions and only finite clusters 
of nonsuperconducting junctions. A more accurate estimate of a would require a precise determination of J c , many 
averages over disorder and long simulations due to the divergent relaxation time near J c as discussed below. Despite 
the uncertainties in the estimate of a, the behavior in Fig. 2b for the largest system size suggests that a could be a 
universal critical exponent independent of the degree of disorder parametrized by p as long as p c < p < 1. 
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FIG. 2. (a) Critical current J c and linear resistance as a function of p for L — 64. (b) Exponent a of the power-law 
behavior E ~ (J- J c ) a . 

At or sufficiently close to p c , where the percolation correlation length is the dominant length scale, one expects a 
different behavior. In fact, by matching the scaling of the IV curves below and above p c , it has been shown thatfl at 
p c the power-law exponent a = 1 + s/[v(d — 1)] with J c = which gives a ~ 1.98 in d = 2. We find a = 2.1(2) for 
L = 64 (not indicated in Figure 2b) which is consistent with this value. 

Some insight into the possible universal behavior for p > p c can be obtained by regarding the onset of resistive 
behavior for J > J c as a dynamical critical phenomenon driven by the external current, where power-law behavior 
appears naturally as a result of scaling. The required scaling assumptions are similar to those proposed by Fisher et 
aZ.tj in a different context. Above the critical current, the superconducting coherence length £ is finite, leading to 
resistive behavior. We assume that it diverges as a power law £ ~ ( J — J C )~ UI near J c , where Vi is a critical exponent 
characterizing the current driven transition. From the definition of the electric field we have E = —dtA, where A is 

the vector potential which enters the Hamiltonian of the JJA in the dimensionless form A ■ dr. The typical time 
scale is given by the relaxation time which diverges near J c asT~^ z , where z is a dynamic critical exponent. From 
dimensional analysis we then expect the scaling E ~ 1/(£t) ~ and a power-law behavior of the current-voltage 

curve 

£~(J-J c ) a , a=[z + l)v I , (5) 

above J c . Sufficiently close to the percolation threshold, there is another characteristic length scale, £ p ~ (p — p c ) ~ Vp , 
the percolation correlation length, and one expects a crossover to another critical behavioifi. We have only considered 
the case of a single relevant length scale and focus in the regime p > p c . 

The exact values of the critical exponents vj and z are not known but the following qualitatively estimate of these 
exponents appear consistent with the numerical results. We expect that at a characteristic current density J > J c , 
phase coherence can change significantly in a correlation volume of order £ d . In this volume, the typical variation of 
the phase is ~ 2-7T. By comparing the coupling term of the external current to the phase gradient J (J — J c ) which 
appears in the continuous version of the JJA Hamiltonian and the quadratic approximation to the Josephson energy 
term K* J (V#) 2 , where K* is an effective stiffness, one finds that (J — J c ) should scale as (J — J c ) ~ l/£ and so vj = 1 
. Thus, the only remaining parameter is the dynamical exponent z. If we neglect non-linearities and assume relaxation 
dynamics as in time-dependent Ginzburg-Landau models, we expect z — 2, and therefore a — (z + = 3. On the 
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other hand, recent workcHl suggests that for the dynamics described by Eq. (||), where there is current conservation 
at each lattice site, the dynamical exponent is z ~ 0.9 and so a = 1.9. The data in Fig. 2b for the largest system size 
L = 64 is intermediate between these two values. Since in Eq. (||) a depends both on vj and z, we need additional 
results to test this analysis. 

We have performed a finite-size scaling analysis at a fixed value of p = 0.7 to verify the scaling behavior of Eq. 
(|J) and extract an independent numerical estimate of vj . In a finite system, the correlation length is limited by the 
system size L and finite-size scaling leads to 

EL a l VI = /((J — J c )L l l Ul ), (6) 

where / is a scaling function. From Eq. (^), all data in the scaling plot EL a l Vl x ( J — J C )L 1 / 1 ' 1 should collapse on to 
the same curve if v and a are chosen correctly as shown in Fig. 3 for system sizes ranging from L = 16 to 64. We find 
that a reasonable scaling behavior is obtained for J c ~ 0.19(2), v\ ~ 1.1(3) and a ~ 2.4(2) where the error estimates 
are obtained by averaging various estimates of J c , vi and a which give equally acceptable scaling plots. Using the 
relation in Eq. (|^), we find z = 1.2(6). The value of vj is in agreement within the errors with the one predicted above 
but z is not accurate enough to allow any comparison. Improved estimates and further analytical work are necessary 
for a detailed study of the critical behavior. The close analogy to other critical phenomena near threshold^ may 
provide an interesting approach and eventually identify the relevant universality class. 
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FIG. 3. Finite-size scaling plot of the power-law behavior E ~ ( J — J c ) a for p = 0.7 at zero temperature, using J c = 0.19, 
Uj — 0.9 and a = 2.5. 

It is interesting to note that if we assume Vj = 1 in Eq. (^) , the result a = z + 1 is similar to the one inferred by 
PresteniH based on an analogy between the onset of dissipation at the critical current and the resistance of a random 
resistor network which leads to a = t + 1, where t isihe conductivity exponent of a mixture of resistors and insulators. 
This would suggest z = t. In two dimensions, whereEZIt ~ 1.3, this gives z = 1.3 and a = 2.3 which is in fact consistent 
with our numerical estimates. 

We now turn to finite temperature effects. |We have only studied the behavior at percolation threshold where the 
superconducting transition is known to occuiliU at T — 0. Again, this is the simplest case where there is a single 
dominant length scale in the system which at finite temperature is the thermal correlation length £j\ This correlation 
length diverges for decreasing temperature as £t oc \/T Ut where vt is the thermal correlation length exponent of the 
diluted XY model which is isomorphic to the J J A at zero current. One can anticipate that .th.r .increasing correlation 
length will have important effects in the nonlinear resistance for decreasing temperaturesoE3. At any T > 0, the 

linear resistance, Rl = Uruj >o E / J , is nonzero since the superconducting correlation length is finite and is expected 

to be thermally activated. However, in presence of a finite current density J, an additional length scale Li oc kT / J is 
set by the external currentllil due to temperature fluctuations. This can be obtained by comparing the extra energy 
arising from the coupling to the external current, JLiAO within a length scale L\ and for a typical phase variation 
A9 ~ 2tt, to the thermal energy kT. For £y < Lj, which holds for sufficiently small J at any finite T, the linear 
resistance is basically unchanged since the smaller length scale £p dominates the activation energy. However, for 
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current densities larger than J n i oc T 1+VT nonlinear behavior sets in as Lj < £r in this case. So, the range in J 
where E/ J is roughly a constant should decrease with temperature and the characteristic current density J n \ where 
it crosses over to nonlinear behavior decreases as a power law T 1+UT with a universal exponent. Associated with the 
divergent correlation length £r one also defines a relaxation time r that owing the zero-temperature transition does 
not follow the usual form r oc £ z and can have an exponential temperature dependence. Since the electric-field scale 
as E ~ 1/(£t), the current density as J ~ fcT/£ and using £ = \jT VT the nonlinear resistance behaves asEJ 

E 1 . J . . . 

in two dimensions, where 17 is a scaling function. If the linear resistance Rl is finite for any T > then 5(0) is a 
constant which can be set to unity, g(0) = 1. When the nonlinear resistance E/J is normalized by its linear value i?z, 
at the same temperature one can then write 

It is clear from Eq. (^|) that the characteristic current at which nonlinear behavior is expected to set in varies as 
jii+i/ x ag men tioned before. 

The nonlinear resistance E/J obtained numerically at p = p c for the largest system size L = 64 is shown in Fig. 4. 
The curves shows the expected behavior. For small current densities J, there is a linear contribution where E/J tends 
to a finite value which depends on the temperature. This is more clear for the highest temperature T = 0.7 where the 
range of J in which E / J is roughly a constant is more pronounced. For increasing current densities it crosses over to 
a nonlinear behavior. As temperature decreases nonlinearity appears at smaller currents and the linear behavior is 
less clear. For the lowest temperature T = 0.3 the linear behavior presumably occurs at current densities smaller than 
J = 0.02 but numerical calculations in this range requires very long equilibration times which prevent us to confirm 
this behavior. In fact, as discussed below, the relaxation time r increases very rapidly (possibly exponentially) with 
decreasing temperatures. 
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FIG. 4. Nonlinear resistance E/J as a function of temperature for L — 64 at p = 0.5. Continuous lines are a guide to the 
eyes. 

The temperature dependence of the linear resistance estimated at the lowest current is indicated in Fig. 5a. It seems 
to be consistent with an activated behavior with an energy barrier ~ 0.91(3). Our data is not accurate enough 
and the temperature range is too limited to exclude more complicated behavior as a temperature dependent Ef,(T). 
The links-nodes picture of the percolation clusterEZI might suggest a single-junction behavior which in fact gives an 
Arrhenius behavior at sufficiently higher temperaturesEJ but this would give a much larger barrier of E^ = 2 and 
requires that the Josephson coupling of the effective single-junction is not renormalized. In addition, the scaling form 
in Eq. (pi) which is found in our case as discussed below, does not hold for the single junction as can be verified from 
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the closed-form solution of the current-voltage relation!!. It is unclear at the moment what is the appropriate model to 
describe the temperature dependence of the R L . In any case, if the apparent exponential behavior of Rl holds down to 
very low temperatures it implies, from Eq. (]?]), that the relaxation time diverges exponentially r ~ exp(Eb/T)/T and 
so the current- voltage characteristics at very low temperatures and currents may be inaccessible by direct simulation. 

In Figure 5b, we show the temperature dependence of J n i for the data in Fig. 4 in a log-log plot. It is consistent 
with the power law behavior J n i oc T 1+VT and provides a direct estimate of the thermal exponent vt — 1-2(2). To 
estimate J n i we defined the crossover to nonlinear behavior as the value of J where E/JR L starts to deviate from 
a fixed value equals to 1.2. The slope in the plot of Fig. 5b should not dependent on this value as long as it is not 
too large compared to unity and we checked that other choices gives the same results within the error estimates. In 
Figure 6, we show the scaling plot of the data in Fig. 4 according to the scaling behavior of Eq. (^). The scaling 
plot is obtained by adjusting a single parameter that gives the best data collapse. This is consistent with the 
scaling behavior discussed above and gives an independent estimate of vt ~ 1-4(2). From these results we obtain a 
final estimate of the thermal critical exponent vt — 1-3(3). This valu&,of vt is in fact consistent with the thermal 
correlation exponent of the diluted XY model at percolation threshold^, v = 0.98 — 1.03. 
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FIG. 5. (a) Arrhenius plot for the temperature dependence of the linear resistance Rl, estimated at J = 0.02, for L = 64 
at p = 0.5. (b) Crossover current density J n i where nonlinear behavior appears as a function of temperature. The slope gives 
and estimate of 1 + vt = 2.2(2) 



Our analysis for the temperature effects on the nonlinear resistance is strictly valid at p c , where the percolation 
correlation length £ p is infinite but the thermal correlation length £t is finite. In order to compare to the available 
experimental data of Harris et aln on artificial arrays for p close to percolation threshold, which is in fact consistent 
with a vanishing transition temperature, we must take into account the competing effects of £t and £ p which is a 
more complicated problem. However, sufficiently close to p c the analysis should still be valid at high temperatures 
when £y << Cp- Unfortunately, the scattering of the data at small currents and the limited range of temperatures 
where linear resistive behavior is apparent prevent us to perform the same scaling analysis as described above. 
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FIG. 6. Scaling plot of the data in Fig. 4 for vt = 1.4. 



IV. CONCLUSION 



In summary, we have studied the IV characteristics in a model of resistively shunted two-dimensional diluted JJA, 
at zero external field, by numerical simulations and scaling arguments. At T = 0, the IV curves show power-law 
behavior above a critical current density which decreases with dilution. The power-law behavior follows from a simple 
scaling analysis which leads to a = (z + where z is the dynamical exponent and Vj is the superconducting 

correlation length exponent. Numerically we find vj = 1.1(3) , consistent with a scaling argument which gives vi 
— 1 , and a = 2.4(2) . The value of a is in agreement with the relation a = t + 1 , in two dimensions, which has 
been proposedtHI in relation to granular high-T c materials in zero field. At the percolation threshold and finite T, 
the results are consistent with the scaling behavior of a T = superconducting transition. Crossover to non-linear 
behavior appears at currents proportional to T + UT , where vt is a correlation length exponent for the diluted XY 
model at percolation threshold. The behavior at percolation threshold is analogous to the zero-temperature vortex 
glass of disordered superconductors in a magnetic field, except for the value of Vt . For experiments in arrays and 
granular high-T c materials, this behavior clearly demonstrates the importance of carefully comparing the expected 
power-law behavior of IV characteristics resulting from field-induced effects to the zero-field case. 
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